Equilibrium roughening transition in a ID modified sine-Gordon model 
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We present a modified version of the one-dimensional sine-Gordon model that exhibits a thermo- 
dynamic, roughening phase transition, in analogy with the 2D usual sine-Gordon model. The model 
is suited to study the crystalline growth over an impenetrable substrate and to describe the wetting 
transition of a liquid that forms layers. We use the transfer integral technique to write down the 
pseudo-Schrodinger equation for the model, which allows to obtain some analytical insight, and to 
compute numerically the free energy from the exact transfer operator. We compare the results with 
Monte Carlo simulations of the model, finding a perfect agreement between both procedures. We 
thus establish that the model shows a phase transition between a low temperature flat phase with 
intriguing non trivial properties and a high temperature rough one. The fact that the model is one 
dimensional and that it has a true phase transition makes it an ideal framework for further studies 
of roughening phase transitions. 



PACS numbers: 81.15.Aa, 68.35.Ct, 68.35. Rh, 05.40.-a 



I. INTRODUCTION 

The two dimensional (2D) ordered sine-Gordon model 
is today a fairly well understood problem (see, e.g., 0,13; 
0, II, la lg|)- However, the random version of the model, 
where quenched disorder is introduced, is far less under- 
stood and still subject of discussion. Since the super- 
roughening transition (see Sec. ^ below for a definition) 
was introduced in 1990 7j, there has been no theoreti- 
cal agreement about its nature and the properties of the 
super-rough phase (see jS] for a review, see 'yi| for more 
references). Large-scale numerical simulations or exact 
optimization results [l3|ll|ll|ll|ll|ll|ll|ll|ll, 
were not enough to solve the question, due to its highly 
demanding computational character. To circumvent this 
problem, in 19] we proposed a modification of the one 
dimensional (ID) model, in order to have a less demand- 
ing problem that could give us information about the 
super-rough phase. In the present work, we proceed to 
a detailed characterization of the ordered version of the 
model in |l3 and the roughening transition it presents. 
Having such a thorough analysis will not only serve as 
grounds for our results on super-roughening [Ift] but will 
also be relevant from a much more general viewpoint, as a 
case study for ID phase transitions and as an alternative 
way to obtain information about complicated problems 
in higher dimensional systems. 

Indeed, the first of the two goals above may seem ques- 
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tionable in view that the subject of ID thermodynamic 
phase transitions, defined as non-analyticities of the free 
energy, has been for a long time excluded from the at- 
tention of the community. This exclusion comes from the 
'public knowledge' that phase transitions can not occur 
in ID systems with short range interactions. However, 
this general belief has risen due to the misunderstanding 
of van Hove's theorem |23, |21| and abuse of Landau's 
22] argument about the entropic contribution of domain 
walls to the free energy. In fact, there are many known 
examples of this kind of transitions (see |23 for a com- 
prehensive study of the matter) , although most of them 
have been hidden using language tricks that have made 
us speak about 1-1-1 dimensions. This has been the case, 
for instance, with a number of models of the so called "2D 
wetting" that are in fact ID in the mathematical sense, 
as J4, _25], just to give a couple of examples. In this 
sense, our work is yet another piece of detailed evidence 
about the existence of ID phase transitions. In addition, 
our model has immediate applications, such as crystalline 
growth over an impenetrable substrate, or "2D wetting" 
favoring the formation of layers of the condensed phase. 
On the other hand, as mentioned in the previous para- 
graph, it is clear that the same approach can be used for 
the study of other 2D problems, that may allow the for- 
mulation of a ID counterpart as we present in this paper 
for the sine-Gordon model. 

With the above objectives in mind, the paper is or- 
ganized as follows: Section ^1 introduces our model 
by discussing in detail its predecessors, the Burkhardt 
and the sine-Gordon ones. Subsequently, in Sec. IIIII 
we present the transfer operator formalism and develop 
it into the pseudo-Schrodinger-equation approximation 



that predicts a phase transition for the model. We thus 
obtain analytical expressions for the magnitudes of inter- 
est at low and high temperatures. Then in Sec. II VI we 
solve numerically the transfer operator problem, show- 
ing the existence of the phase transition and computing 
thermodynamical magnitudes such as the specific heat. 
In Sec. we present the results of Monte Carlo parallel 
tempering simulations of the model, compare them with 
the preceding results and discuss the non-trivial behav- 
ior in the flat phase of the model. Finally, in Sec. IVll 
we discuss the consequences of all this as well as further 
implications of our results. 



II. MODEL DEFINITION 

In order to properly introduce and motivate our model, 
we find it convenient to review in some detail the two 
previously proposed ones on which it is based, namely 
the model introduced by Burkhardt in 1981 |2y| in the 
context of wetting, and the sine-Gordon model, a widely 
applicable model representative of a variety of physical 
systems (see Q and references therein). Beginning with 
the first one, the Hamiltonian of Burkhardt's model is 
given by |26l |: 



'H = 



N 



-h,\- 



u{h,)y 



(1) 



and defines a continuous counterpart of the models pro- 
posed W Chui and Weeks [23 and van Leeuwen and Hil- 
horst J2gj in 1981. We are interested in the version of 
the model with positive values of the variables {hi > 0). 
U{hi) is a function with a positive constant value Uq 
for hi < R and zero for h > R. The variables hi can 
be seen as heights over a substrate (located at /i = 0), 
defining all together an interface. This model is exactly 
solvable because its transfer integral equation can be ex- 
actly mapped |26| to a Schrodinger equation, formally 
a quantum square well problem in ID with the square 
well potential at the edge of the system. From quantum 
mechanics |29|| we know that this potential has a bound 
state solution for a well deep enough. In Burkhardt's sta- 
tistical mechanical problem, the depth of the well of the 
resulting Schrodinger equation depends on [3, the inverse 
temperature. Hence, for low enough temperatures, the 
quantum bound state maps to an interface trapped by 
the potential, and therefore the interface is flat, in the 
sense that its width is finite. Above the critical temper- 
ature of the model, the bound state disappears and the 
interface dcpins from the potential and its width diverges: 
it becomes rough. 

The other pillar on which our model stands is the ID 
sine-Gordon model, whose Hamiltonian is 

N J 

n = J2 {2^^^+^ ~ ''')' + M^- cos(/iO]}, (2) 



where now the values of the variables are unrestricted 
(—(X) < hi < cxd). Its 2D version is very interesting 
because it can describe a number of different physical 
problems 9] and because it presents a roughening phase 
transition. Again in the language of hi being the height 
of a surface, this transition takes place between a high 
temperature rough and a low temperature flat phase (we 
will define more precisely below what we understand by 
rough and flat). In the rough phase, the roughness (also 
to be defined below, but in the surface language can be 
thought of as the surface width) of the system scales as 
InL, the logarithm of the system size. The roughening 
transition is modified by the addition of disorder to the 
system: when the cosine potential is changed adding a 
quenched disorder h^, making it Vb[l — cos{hi — hi)], a 
super-roughening transition arises, characterized by the 
fact that the low temperature phase is no longer flat. 
The super-roughening transition and specially the low 
temperature phase are poorly understood. One of the 
features that seem to be accepted about this super-rough 
phase is that in it the roughness scales as In^ L, so in this 
sense it would be even rougher than the high tempera- 
ture rough phase (hence the name super-rough). Unfor- 
tunately, the ID version of the sine-Gordon model, much 
easier to study analytically and numerically, is of no help 
to shed light on this problem, as long as it has been rig- 
orously proven 30] that it can not have a true thermo- 
dynamic transition (although it docs have an apparent 
one for any finite size systems, even extremely large ones 
P). 

In order to better understand the phenomenology of 
the 2D version of the model, in [l^ we introduced a new 
model containing the features of Burkhardt's and sine- 
Gordon models, in order to retain the most interesting 
characteristics of both of them: the phase transition of 
the first one and the periodic potential of the latter. The 
rationale for this approach was that if we had a ID model 
with such a phase transition, we could consider its dis- 
ordered version and check whether or not it reproduces 
the features of 2D super-roughening. We indeed carried 
out this program in 19], but a key question remained, 
namely whether or not the basic, ordered, ID model had 
a true thermodynamic phase transition or not. Only if 
the answer of this question is positive will the approach 
in 19] make sense. Although our model is specifically 
designed to exhibit this phase transition, and hence the 
transition itself would not be surprising, we must prove 
hat the model behaves as expected: the fact that the 
model ingredients suggest that it will indeed have a tran- 
sition by no means warrants its existence. In addition, 
the main purpose of this paper is the characterization of 
the low temperature phase, which will show novel non 
trivial behavior as we will see below. 

The Hamiltonian for our model, which we called the 
Burkhardt-sine-Gordon model, (BsGM hereafter), is: 



n = 



N 



- hi 



Vr 



o-V{h,)Y 



(3) 



where 



III. ANALYTICAL RESULTS. 



V{x) 



[1- 

CX) 



{x)]-^U{x) ifa;>0, 



if a; < 0. 



(4) 



The choice for a quadratic coupUng instead of an ab- 
solute value one as in [2g is to make our model as close 
as possible to the original sine-Gordon model. In addi- 
tion, the gaussian fluctuations of the heights that this 
coupling implies can be simulated with higher efficiency 
using a heat bath algorithm [S^, 123 ■ We impose no ex- 
plicit restriction over the values of hi, it is the value of 
the potential for hi < what forces the variable to take 
only positive values. This unlimited range of definition 
of the variable will be useful for the formal operations we 
will perform. We see now that U{x) can be seen as an at- 
tractive potential binding the interface to the substrate. 
The cosine potential will favor the growth forming lay- 
ers at a distance 2tt of each other. For definiteness we 
choose the parameters of the model to be Vb = 1, C/q = 2 
and R — 2tt. In that way, the substrate will attract the 
two first layers. We have also performed simulations with 
different values of the parameters. In that way, we can 
change the number of attracting layers, or the critical 
temperature, but no new qualitative behavior is found. 

In view of the above considerations, we expect the 
BsGM to have a phase transition between a flat (or 
pinned) phase at low temperatures and a rough (or de- 
pinned) one at high temperatures. That is exactly what 
we needed in order to compare to the results on 2D super- 
roughening in disordered sine- Gordon models 19] . How- 
ever, in that previous work, we did little more than pro- 
viding plausibility arguments and simulation evidence for 
the existence of such a transition, hence the necessity of 
the detailed, much more rigorous work presented here. 
To characterize the transition, we define the roughness 
or interface width, w, as: 



where 



-Uh^^-w 



h = ±yh 



(5) 



(6) 



is the mean height, and averages (• • • ) are to be under- 
stood with respect to a statistical weight given by the 
Gibbs factor, e~'^''^ , at equilibrium at a temperature T. 
Then we say that an interface is flat when w is finite and 
does not depend on the system size, N. In the rough 
phase, the interface width grows with N and diverges in 
the thermodynamic limit, w —^ oo as N —^ oo. Addition- 
ally, in the remainder of the paper we will look at other 
possible indicators of the transition, such as the free en- 
ergy, the specific heat or the full correlation function. 



A. Transfer integral technique. 

The following discussion of the transfer integral (TI) 
technique follows that in 33] for the sine-Gordon model. 
The classical canonical partition function of the BsGM 
[Eq. ©] can be written as: 

/CXD /"OO /.OO 

d/ii / dh2--- dh^e-^^, (7) 

-OO J — OO J — OO 

(3 being the inverse temperature in units of the Boltz- 
mann constant. Note that we could have written the 
integrals in the range [0,cxd), but our definition of V{hi) 
makes this unnecessary. In what follows, periodic bound- 
ary conditions 

hi = hN+i (8) 

are assumed, so that Eq. lO can be replaced by 

/OO /-OO 

d/li---/ dhM+l(i-^^5[hi-hM+l). 
-QQ J —OQ 

(9) 
To evaluate Zn{/3) we proceed as follows. The S func- 
tion is represented as an expansion in a set of complete 
orthonormal functions ipn{h): 



5{h-h')^Y.^l{h)^,,{h'). 



(10) 



The functions (^„ are chosen to satisfy the the TI equa- 
tion: 

dheM-PyQK{h,h')]^,,{h) = exp(-/31/oen)'^«(/i') 

(11) 
where 

K{h, h') ^~{h- h'f + \[V{h) + V{h% (12) 

Lpn is an eigenfunction of the TI equation with associated 
eigenvalue exp(— /3Vben), and 



r(/3) = exp[-/3VoX(/i, h')] 



(13) 



is the transfer operator of the model. Using this we can 
rewrite the partition function: 

/OO 
dh^*^{h)^n{h) = 

^Y^eM-PVo^uN). (14) 

n 

In the last step we have used the orthonormality of the 
ipn- The orthogonality and completeness of the eigen- 
functions are guaranteed by the Sturm-Liouville theory 
for Fredholm integral equations with a symmetric kernel 
|34l |. for which Eq. H12|l is an example. 



Since the single site potential [Eq. |0J] is bounded from 
below, the eigenspectruni is also bounded from below, 
and we denote the lowest eigenvalue by eg. This cor- 
responds to exp(— /3eo), the maximum eigenvalue of the 
transfer operator H13|l . In the thermodynamic limit, the 
free energy per particle is then given by 

/ = -ksT lim 1 liiZNil3) - yoeo- (15) 

N^OQ iV 

From this result, other thermodynamic properties can 
now be derived, i.e., internal energy per particle 

e=^(H) = /-T|^ (16) 

and specific heat at constant volume (length) 
de d^f 

In [221 it is also shown that in the thermodynamic limit 
canonical averages are given by the expression 



where 



{9{h^)) = 



\Mh)\^9ih)dh, 



(18) 



what means that |(/3o(^)P can be understood as the prob- 
ability density for the variables hi. 

We have been unable to exactly evaluate the free en- 
ergy ^^ for the BsGM [Eq. ©]. As an alternative, in 
this work we will proceed in two different ways: through 
the pseudo-Schrodinger-equation approximation associ- 
ated to the TI equation and computing numerically the 
eigenvalues of the transfer operator. In what follows we 
discuss the former approach, as well as the approximate 
analytical results that can be obtained from it. The nu- 
merical study of the full transfer operator deserves a sep- 
arate treatment and is reported in Sec. IIVI 



B. Pseudo-Schrodinger-equation. 

Defining 

Mh) = eM-PVo^V{h))^„{h) 



(19) 



and using the identity 
1 



dy exp 



2t 



{x~y)']]f{y)^ 



the TI equation (|ll|l may be rewritten in the form 



(20) 



exp 



exp 




Mh) (21) 



V,, 



1 , (131 
2VqI3 ^\2tt 



(22) 



In the limit of strong coupling ( J/Vq — > oo as we keep VqJ 
constant, see [33] for details; it can also be understood 
as a continuum limit if we include the lattice constant as 
an explicit parameter of the model) this equation can be 
expanded to obtain, to first order in Vq/ J: 



1 



2p^VoJdh^ 



V{h) 



Mh) = {er,-Ve,)Mh). (23) 



This is the Schrodinger equation for a square well with 
a superimposed cosine potential. The square well alone 
is exactly solvable, and for low values of Vb the cosine 
term can be treated with perturbation theory. Eq. H23|l 
is already an approximation for our model; from Eq. 
(|21|l we can see that it is valid when y/Vo/J <^ (3 and 
\/VqI J ^ 1//3 . If we take Boltzmann's constant as 
unity, the predictions of this equation are expected to 
hold quantitatively only in the the temperature region 
^JVq/J < T < 1/^JVq/J. For this interval to make 
sense, ^JVqJJ has to be small enough. For instance, for 
the values of the parameters we use (J = 1, Vq — l)j 
there is no temperature where the pseudo-Schrodinger 
equation is quantitatively accurate. However, the quali- 
tative picture this equation yields is completely valid and 
describes correctly the phenomenology of the model. In 
the quantum mechanical problem, for some values of the 
parameters of the model we have a bound state, that 
disappears as we change the parameters. In our statis- 
tical mechanical problem, fixing all the parameters ex- 
cept the temperature will give us a thermodynamical 
phase transition between a flat phase at low tempera- 
tures, pinned by the square well potential, and a rough 
phase at high temperatures, where the interface has de- 
tached itself from the substrate's attraction. This is the 
same scenario Burkhardt found in j2Q| ; the change of the 
absolute value coupling for the quadratic one and the 
addition of the cosine potential modify the quantitative 
aspects of the phase transition, but not the qualitative 
ones. Of course, these new features in our model will 
give rise to new phenomena in the flat phase's behav- 
ior. Anyway, if we make a further rough approximation 
and dismiss the sinusoidal part of the potential in Eq. 
l|^ . we are left with exactly the Schrodinger equation 
of a semi-infinite square well. From elementary quan- 
tum mechanics ^3 (see also 0| for the application to 
Burkhardt 's model), we know that the spectrum of this 
equation presents a continuum of scattering states. For 
appropriate values of the parameters (that in the statis- 
tical mechanical problem means T < T^) there are one 
or more bound states. As T -^ T^ ^ the gap between the 
strongest bound state and the first scattering state varies 
as: 



Ae oc (T, - Tf 



(24) 



V(h) 




we have for the internal energy: 



FIG. 1: Approximation of the potential inside the square well 
by a (f>'^ potential. The continuous line is the BsGM potential 
between and 27r. The dashed line is the (f>'^ potential we use 
to approximate it. 



The quadratic temperature dependence of the gap in 
Eq. (|24|l is responsible for the finite jump in the specific 
heat of the model. We will find this in the computation 
of the specific heat both from the numerical transfer op- 
erator and from Monte Carlo simulations. In Figure|31we 
show the gap between the two first eigenvalues computed 
from the exact numerical transfer operator; the quadratic 
behavior predicted in Eq. (|24|l is evident as T — > Tc~ . 

For the rest of this work, without loss of generality, 
we will take the coupling constant J = 1. We can do 
this because the effect of changing J can be taken into 
account rescaling Vb, Uq and the temperature (and also 
the time scale, but in this work we will deal only with 
equilibrium properties). 



C. Low and high temperature approximations 

For low enough temperatures, it is a good approxima- 
tion to suppose that all the heights fall inside the square 
well potential. For a value of the width of the well of 
R = 2iT, inside the well there exist two minima of the 
cosine potential. In that case it is reasonable to approx- 
imate the potential by a t/)^ one, see Fig. Q] The good 
features of this choice are that the </>'' potential repro- 
duces the two potential minima and that it bounds the 
system to them, as it grows to oo as /i — > ±cx). Note 
that if we restrict ourselves to only one minima of the 
cosine potential, a parabolic potential will be enough to 
reproduce the leading term. To mimic the potential in 
our problem, this (p'^ potential has the form (for Vq — 1 
and [/o = 2): 



V(/i) 



{h 



(h- 



47r2 



T 



(25) 



In |33| we find values for some thermodynamic properties 
of a low temperature expansion of the cj)'^ model. Thus, 



T 
2 



36T^ 



15-237r2' 



and for the specific heat: 



cv 



72T 



15 • 237r2 



(26) 



(27) 



We will see that, at low temperatures, the system chooses 
to be in one single potential minimum of the two dis- 
played in Fig. ^ In fact, this assumption is implicit in 
the calculation that lead to Eqs. l(^ and (|77|l (see JSJ for 
details). This calculation approximates the (p'^ potential 
by a parabolic one {V{h) — Voh'^), and then introduces 
the higher order corrections. 

The same procedure can be used with the sine-Gordon 
model instead of the 0^ model. It also seems a reason- 
able choice to approximate the T — > regime using this 
potential. In the end, as both models have the same lead- 
ing term, the differences between them will be small. We 
will compare the expressions arising from both of them 
with the results of our simulations, and find that both of 
them describe remarkably well physical magnitudes when 
T —f 0. From 33], we have the following expressions for 
the low temperature sine-Gordon model: 



e^-+2 



T 



cv = - + 2 



2T 



T 



3T2 



(28) 



(29) 



Both these approximations suppose the system is 
trapped in a single well of the potential, and it can be 
seen that this implies that the system is in a fiat phase 
|31| . So agreement with these results is a signal of a flat 
phase. 

Restricting ourselves to the lowest order approxima- 
tion for vanishing temperatures, that is, a flat system 
trapped in a single parabolic potential, it is straight- 
forward to calculate the roughness and correlation func- 
tions, as was done in |3lj| . The parameter of the parabolic 
potential has to be Vb -I- C/q . For the roughness we obtain: 



\T)^ 



T 



V(2 + Vb + [/o)2 - 4 
We define the height-difference correlation function as 



(30) 



Cir)^(^Y.^h,-h,+r] 



(31) 



It can be shown that the parabolic potential approxima- 
tion yields for it: 



C(r) = 



2T 



y/{2 + Vo + U^Y - 4 



(1 - C,(r)) (32) 



where 



Ccir) 




Vo + Uo 



1- Wl 



2 + Vo + Uo 



(33) 

In the asymptotic hmit r ^ oo, Cc(?') —> 0, and we have 
that C(oo) = 2w^. 

In the high temperature phase, the potential effectively 
vanishes and we are left with the quadratic coupling 
alone: this is the Edwards- Wilkinson model [33 • The 
predictions for the internal energy (e = T/2 + constant) 
and the specific heat (cy = 1/2) are expected to hold in 
the rough phase of our model. However, the prediction 
for the interface width is not so accurate: the existence 
in our model of an impenetrable substrate changes the 
statistics of the rough interface. 



IV. NUMERICAL TRANSFER OPERATOR 
RESULTS 
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■ ■ ■ Third eigenvalue (exp(-Pe,)) 




FIG. 2: Three first eigenvalues for M = 4096 and Ah = 1/32. 
Inset: Difference Ae = exp(— /3eo) — exp(— /3ei) vs. T. The 
minimum gives the temperature of the phase transition. 



The eigenvalue problem in Eq. (|ll|l can be solved dis- 
cretizing the transfer operator in Eq. (|13|l and evaluat- 
ing numerically the eigenvalues of the resulting matrix 
(see [33, ISa 123; see |3g| for a detailed account). The 
relevant parameters of the discretization of the opera- 
tor are Ah, the discretization length, and M , the num- 
ber of points considered, that is, the size of the matrix. 
From them we obtain immediately the interval where 
the discretized variable takes values, [0,hmax], where 
hmax = {M — 1) • Ah. The two sources of error of this nu- 
merical procedure are the discretization of the real vari- 
able h and the cutoff of the variable range at hmax- In 
the limit Ah — > and MAh — > 00 (that is, hmax —^ 00), 
this numerical approach is exact. 

A thermodynamic phase transition takes place when 
there is a non-analyticity in the free energy. We have 
seen in Eq. H15|) that in the thermodynamic limit the 
free energy is determined by the largest eigenvalue of the 
transfer matrix. As discussed below Eq. (|24|l . the vanish- 
ing of the gap between the largest two eigenvalues leads 
to a singularity. To find the point of a phase transition, 
we have to find a minimum of the gap and show that the 
minimum goes to zero as we increase M. 

In Figure HI we show the first three eigenvalues of the 
discretized transfer operator with our standard set of pa- 
rameters, Vq = 1, t/o = 2 and R = 27r. We clearly 
see that the first two eigenvalues become very close near 
T w 10. In the inset we show the minimum of Ae that in- 
dicates the temperature of the candidate transition. The 
slope of eg does not change discontinuously at T„n the 
temperature of the minimum, so the transition will be 
continuous and not first order. In Figure Owe show the 
gap between the two first eigenvalues for a range of ma- 
trix sizes, keeping Ah fixed. We see that as M increases, 
the minimum value of the gap becomes closer to zero. 
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FIG. 3: Ae for different matrix sizes as indicated in the plot. 
The discretization is Ah = 1/8. Inset: the same figure with 
the Ae axis in logarithmic scale. We see that as M becomes 
greater, Ae goes quadratically to its minimum as T -^ T^ , 
as shown for M = 1536 using a quadratic fit. This is exactly 
the prediction of Eq. 11241 . 



In Figures 01 and we perform a finite-size scaling to 
check that the minimum of the gap, Aemin, goes to zero, 
and how the different temperatures for the minimum go 
to the critical temperature, Tc- We see, as observed 
in pjq for a different model, that both Aemin a-nd T^ 
scale with M~^ when we change M keeping Ah fixed. 
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FIG. 4: Minimum value of the gap for different discretization 
values and matrix sizes, as indicated in the plot. 
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FIG. 5: Critical temperature for different discretization values 
and matrix sizes as indicated in the plot. 



Of course, this scaling is supposed to improve for greater 
matrix sizes, and this aspect is important specially for 
small Ah. In Figure 0] we see how as M increases Aemin 
goes to zero. It may seem contradictory that as we take 
a better (smaller) A/i, the convergence to zero becomes 
worse. The explanation comes from the fact that, as we 




FIG. 6; Specific heat as a function of temperature obtained 
from the discretized transfer operator for Ah = 1/32 and 
M = 4096. 



use a smaller Ah, we need a bigger M to get a correct 
scaling. However, memory limitations of our computers 
sets a limit to the values of M we can use: we can not go 
much further of M — 4096 in a reasonable time. So, to 
get a better estimation of Aemim we use only the points 
with the best scaling. That is what we do for Ah = 1/64, 
where using only the two points of greater M we see that 
the asymptotic value is corrected in one order of magni- 
tude. We can then safely expect Aemin ^ as M~^ — > 
and Ah — > 0. This means that in fact we have a true ther- 
modynamic phase transition, as predicted by the pseudo- 
Schrodinger approximation. The critical temperature Tc 
can be inferred from the data in Figure|Sl The data com- 
ing from the smallest values of Ah are supposed to be the 
best one, and again we have used only the last two values 
for Ah — 1/64 to correct the effects of lack of scaling for 
low M. With the data in the figure, we can estimate the 
critical temperature as Tc = 10.3 in our units. 

We have also computed, using Eq. ((T7|l . the specific 
heat from the numerically obtained eigenvalue. This is 
shown in Fig.|H| The jump of the specific heat at T « 10.3 
is the jump associated with the phase transition. The 
peak at T « 1.4 is a well-known Schottky anomaly (see, 
e.g., (3 and references therein) related to the fact that the 
heights pass from being mostly in one well of the cosine 
potential to expand to different wells. There is an extra 
feature, namely the narrow peak at T « 0.4. If we look at 
the gap Ae between the two first eigenvalues, it effectively 
has a minimum at that temperature, what would make 
us think of an additional phase transition. Furthermore, 
that transition would have a physical interpretation. In 
Fig. we represent the square value of the first eigenvalue 
of the transfer operator, that as we saw in Eq. (|18|1 has 
the interpretation of the probability density of hi. In the 
figure we see that at the temperature of the transition 
(T w 0.77 for the parameters of the figure) the heights 
pass from being almost all in the lowest h well of the 
cosine potential (the potential well with the minimum 
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FIG. 7: Probability density of h for different temperatures 
around tlie narrow peak of the specific fieat for M — 1440 
and hmax = 100 {Ah = 5/72). 



at /i = 0) to being in the highest h one (the well with 
niinimum at h = 27r). This "transition", however, does 
not survive a finite-size study: as, keeping hmax fixed, 
we make Ah smaller, the temperature of the transition 
goes to zero, showing us that it is nothing but a result of 
the discretization and the numerical technique employed. 
Our Monte Carlo simulations will confirm this, as they 
show all the way down to the lowest temperature we have 
simulated, below T = 0.1, that the well preferred by 
the heights is the highest h one (see Sec. |V] and Fig. [T^ 
below). Upon this observation, one question immediately 
arises: if both the first and the second well of the cosine 
potential are energetically equally favorable, why does 
the system choose as the equilibrium one the second? 
The reason is that entropically they are not the same, 
and the configuration of the heights in the highest h well 
has greater entropy. The reason for this is that the only 
escape a height hi has from the lowest h well is going to 
the highest h one (at low enough temperatures at which 
big h differences are very unlikely). But from the highest 
h well it can escape to the lowest h one, or to the next 
cosine well outside the Burkhardt-like square well. So the 
two wells are not symmetrical, and the configurations in 
the highest h one have higher entropy. In that way, what 
we see in the lowest temperature curves in Fig. would 
be in fact a metastable state with higher free energy than 
the true equilibrium one, the heights in the highest h well. 



V. MONTE CARLO SIMULATIONS 

To confirm the conclusions drawn from the analytical 
simulations on the existence of a phase transition we have 
resorted to parallel tempering Monte Carlo simulations 
|3ll Isa . l40| . Representative configurations at a given 
temperature are generated with a heat bath algorithm 
|3lL l3^ , in which new values /i,^ for the height at site i 
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FIG. 8: Internal energy per particle obtained from Monte 
Carlo simulations. Inset: display of the low temperature re- 
gion and comparison with the predictions of Eqs. H26|l and 
12811 . Note that the zero temperature energy is shifted by 
—2 with respect to Eqs. 1261 and 1281 to take into account 
the square well potential. Lines are as indicated in the plot. 
At this scale, the predictions of the sine-Gordon and the ^* 
models (in the inset) are indistinguishable. 



are proposed according to the rule 



/i' 



h.-i-\ + /iM 




(34) 



^ being a Gaussian random variable of zero mean and 
unit variance, and are accepted with a probability 
min[l,e-'5^/'^] with bU = ^{h!^) - V(h,)\. The reason 
to accept or reject using only the potential term in the 
Hamiltonian is that the proposal in Eq. 1341 exactly re- 
produces the quadratic coupling fluctuations, which are 
gaussian. Since that term is already fully included in the 
proposal, we do not need it in the acceptance rate. 

The parallel tempering algorithm then considers simul- 
taneous copies of the system at different temperatures, 
allowing exchange of configurations between them. This 
is particularly efRcient for low temperature configura- 
tions, which are most susceptible to being trapped in 
metastable regions. The simulation starts using a sin- 
gle system copy (replica) at the highest temperature of 
interest. After simulating it we get the temperature for 
the next replica from the energy fluctuations. We repeat 
the same process until we have a set of temperatures 
that covers the whole range of interest. Then we run a 
parallel tempering simulation of all replicas and from it 
get improved values of the temperature set. This auto- 
tuning process continues until we have an almost perfect 
measure of the specific heat, which shows that we are us- 
ing a near to optimal temperature set, and at the same 
time that the different replicas are properly equilibrated. 
After allowing this last temperature set replicas run for 
further equilibration, we start the measuring run. 

The parameters we have used for our simulations, as 
already said, are Vb = 1, f/o = 2 and i? = 27r. We have 
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FIG. 9: Specific heat from Monte Cario simulations, com- 
parison is made witli tlie numerical transfer operator result. 
Error bars of the simulations are of the size of the symbols or 
smaller. Symbols and lines are as indicated in the plot. 
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FIG. 10: Specific heat from Monte Carlo simulations at low 
temperatures compared with the predictions of Eqs. H27|l and 
1291 1. The symbols are simulation results for different system 
sizes as indicated in the plot. Error bars are of the size of the 
symbols. The dashed-double dotted line is the prediction of 
Eq. (I27II and the dotted one the prediction of Eq. (I29II . 



also ran simulations with different values of the param- 
eters without finding qualitative differences. We have 
made simulations for system sizes of A^ = 500, N = 1000 
and N — 2000, although for simplicity we do not present 
results for N = 500. In Fig.|Slwe plot the internal energy 
per particle. We see that the results for both system sizes 
agree perfectly, and that the agreement with the theoret- 
ical predictions for low temperature [Eqs. H26() and (|28|l ] 
is quite remarkable. At high temperature it has the pre- 
dicted slope 0.5, and we see at T ~ 10 the change in the 
slope indicating the temperature of the phase transition. 
Fig. |51 shows the specific heat obtained from the simula- 
tions. We sec that the coincidence between both system 
sizes and the numerical transfer operator result is perfect, 
except in the low temperature region where we have seen 





' 1 ' 1 ' 1 ' 1 


'I'M- 


2000 






/ 






— N=2000 
-- N=1000 


1500 


- 




1000 


- 




500 



■Ill—*' 


¥ - 




02468 10 12 14 01234567 



FIG. 11: Left: squared roughness w^ vs. T. Right: zoom of 
a lower temperature region. Note the perfect overlap of the 
results for the two different system sizes below the transition 
temperature. Inset: yet another zoom of an even lower tem- 
perature region, where we can see the comparison between 
simulation results and the prediction of Eq. 13UII . 
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FIG. 12: Typical interface configuration at low temperatures. 
This one is for the A'' — 2000 Monte Carlo simulation at 
r = 0.0981. 



that the numerical transfer operator introduces the spu- 
rious transition, and in the region of the phase transition, 
where small discrepancies due to finite size effects arise. 
As should be expected, the transition is more abrupt for 
the largest system size, N — 2000. This agreement be- 
tween the results of two completely different approaches 
as the numerical transfer operator and the Monte Carlo 
simulations provides firm grounds to our claims. In Fig. 
llOl we see how the specific heat has an asymptotic behav- 
ior as T — > in agreement with approximations (|27|l and 
(HI. 

As the most important verification of the transition, 
Fig. 111! shows the squared roughness. For temperatures 
above the phase transition, w'^ becomes dependent of the 
system size and diverges with N, showing us that we are 
in a rough phase. Below Tc, the results for both system 
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FIG. 13: Height difFerence correlation functions scaled by the 
temperature from the A'^ = 2000 simulation. Temperatures 
are (from up to down of the greatest value: T =14.0, 10.26, 
9.53, 8.56, 7.80, 6.90, 1.62, 3.99, 1.12, 0.995, 0.836, 0.697, 
0.0981. 



sizes are the same, and as T — > we see the behavior 
predicted in Eq. (|30|l . The step in the roughness between 
T ~ 1 and T ~ 1.5 is an effect of tlic Schottky anomaly 
|8l] | we have already mentioned. Between T ~ 2 and 
T ~ 4, we see a little plateau in the roughness curve. This 
plateau is caused by the dominating part that the kinks 
formed between the lowest h well and the highest h one 
play at these temperatures, while the relaxation of the 
heights in each well as temperature goes down is almost 
screened by the effect of the kinks in the roughness. At 
the lowest temperatures, below T ~ 1, all effects of kinks 
disappear, and the interface is trapped in the highest h 
well in the square potential, as we already noted above 
and show in Fig.^] This is related to the apparent phase 
transition studied in |8]] |. 

In Fig. El we see the height-difference correlation func- 
tion, scaled by temperature, from the simulations with 
N = 2000. All the curves corresponding to temperatures 
higher than Tc collapse to a single curve. This is the ex- 
pected behavior for the high temperature rough phase, 
as the potential term in the Hamiltonian is expected to 
be irrelevant at these temperatures, leaving us only with 
the quadratic coupling, which is the Edwards- Wilkinson 
model [33 that predicts exactly this independence of T 
for C{r)/T, see also 0. The first curve below this col- 
lapse is the curve for T = 10.26. So, from our simula- 
tions we obtain Tc = 10.26, in excellent agreement with 
the numerical transfer operator result. For N = 1000 
(not shown) we obtain Tc = 10.31 and the same behavior 
depicted in Fig. 1131 Finally, note the agreement between 
the prediction of Eq. H32I) and the actual low temper- 
ature correlation functions we find in simulations. We 
see again in Fig. E] the effect of kinks that appeared 
in the roughness between T ~ 2 and T ~ 4: the tem- 
perature scaled height-difference correlation function has 
a non monotonous behavior with temperature between 



T = 3.99 and T = 1.62. In this range the different func- 
tions (without scaling) are almost independent of tem- 
perature, so the scaled functions have higher values as 
we reduce temperature. Note that this behavior only 
appears above certain length scale. At very short scales, 
the effect of kinks has little importance (as we need a cer- 
tain system size to have probabilities of kinks to appear) 
and the relaxation of heights continues with decreasing 
temperature. 



VI. CONCLUSIONS 

We have studied in detail a model first proposed by us 
Il9|. which combines the model proposed by Burkhardt 
in [23 and the well known sine-Gordon model. We show 
here by analytic approximations and by two different nu- 
merical methods (transfer operator and Monte Carlo sim- 
ulation) that it has a continuous phase transition between 
a high temperature rough phase and a low temperature 
flat one. We have characterized the thermodynamics of 
the model, establishing its non trivial behavior in the flat 
phase due to interaction of the two kinds of forces (pe- 
riodic potential and substrate attraction) present in it. 
This gives rise to the existence of a temperature region 
(between T ~ 1.6 and T ~ 4.0) where physical magni- 
tudes of the interface as roughness and spatial correla- 
tions are quite independent of the temperature. In ad- 
dition, our work also stands as a careful study of a ID 
thermodynamical phase transition. While we hope our 
results will stimulate further studies in this field, misun- 
derstood for a long time, we want to add a few caveats 
about how numerical results can lead to misleading con- 
clusions. First, we have seen that the numerical analy- 
sis of the transfer operator produced an artifact which 
looked like a second phase transition in the low temper- 
ature regime. Second, we have shown in a previous pa- 
per f31t| that simulations can yield results reminiscent 
of a true phase transition even for extremely large sys- 
tem sizes, whereas it is rigorously known 30j that such 
transition is impossible. Therefore, it must be born in 
mind that only a judicious combination of theoretical re- 
sults, numerical analysis and simulations may provide 
firm grounds to claims of existence of phase transitions 
in models that are not exactly solvable. This is even more 
important in the case of ID systems, where the debate 
is contaminated by the false prejudices against their own 
existence |23 |. 

Finally, we want to stress that the results we have ob- 
tained on this model suggest a more amenable analytical 
and computationally way to study the properties of mod- 
ified versions of the 2D sine-Gordon model as we did in 
[1^] for the random substrate version. As our model has 
a transition between a low temperature flat phase and a 
high temperature rough one, just like the 2D sine-Gordon 
model without disorder, in that work we showed how the 
addition of disorder to our model can give us insight of 
what happens in the low temperature phase of the 2D 
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random sine-Gordon model. We believe that the same 
ID approach to 2D problems will prove fruitful in many 
other contexts. Its two main advantages are that usually 
ID models are more amenable to analytical treatment 
than 2D ones, and that simulating a ID model requires 
much less computational effort. We hope that many new 
insights will be obtained following this line of research. 
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